Normal Start Functions

Load isofiles

isofiles <-
  file.path(
    "Results",
    c(
      "180305_alkanes_C_GB",
      "180306_alkanes_C_GB",
      "180307_alkanes_C_GB",
      "180309_alkanes_C_GB"
    )
  ) %>%
  iso_read_continuous_flow()
## # A tibble: 6 x 4
##   file_id                  type  func           details                   
##   <chr>                    <chr> <chr>          <chr>                     
## 1 585__Linearity CO2__.dxf error extract_dxf_r… cannot identify measured …
## 2 585__Linearity CO2__.dxf error extract_isoda… `peak`, `start`, `rt`, `e…
## 3 586__Linearity CO2__.dxf error extract_dxf_r… cannot identify measured …
## 4 586__Linearity CO2__.dxf error extract_isoda… `peak`, `start`, `rt`, `e…
## 5 587__Linearity CO2__.dxf error extract_dxf_r… cannot identify measured …
## 6 587__Linearity CO2__.dxf error extract_isoda… `peak`, `start`, `rt`, `e…

Data

File Info

# all file info
isofiles %>% iso_get_file_info()
## Info: aggregating file info from 71 data file(s)

Vendor data table

isofiles %>% iso_get_vendor_data_table()
## Info: aggregating vendor data table without units from 71 data file(s)

Chromatograms

isofiles %>% 
  # fetch a few of the files of interst
  iso_filter_files(parse_number(Analysis) %in% c(616, 632)) %>% 
  # plot just mass 28
  iso_plot_continuous_flow_data(data = c(46)) %>% 
  # make interactive
  ggplotly(dynamicTicks = TRUE)
## Warning: You need the dev version of ggplot2 to use `dynamicTicks`

Analysis

Select relevant data

data_table <-
  isofiles %>%
  # filter the files you want to use 
  # --> exclude CO2 zeros
  iso_filter_files(!str_detect(`Identifier 1`, "CO2 zero")) %>% 
  # --> select only good analyses
  iso_filter_files(parse_number(Analysis) > 572) %>% 
  # get the vendor data table and file info
  iso_get_vendor_data_table(
    select = c(
      # peak info
      Nr., Start, Rt, End,
      # amplitudes and intensities
      `Ampl 44`, `Ampl 46`, `Intensity All`, `Intensity 44`,
      # ratios and deltas
      `R 46CO2/44CO2`, `d 13C/12C`
    ),
    include_file_info = c(
      file_datetime, Analysis, `Identifier 1`, `Identifier 2`, `GC Method`, `Seed Oxidation`
    )
  ) %>% 
  # rename some columns to be easier to work with
  rename(
    Ampl44=`Ampl 44`, Ampl46=`Ampl 46`, # amplitudes
    #Area28=`rIntensity 28`, Area29=`rIntensity 29`, # areas
    Intensity=`Intensity All`, #peak intensities
    R46C = `R 46CO2/44CO2`, d13C = `d 13C/12C`, # ratio and delta values
    rt = `Rt`,
    rt_start = `Start`,
    rt_end = `End`
    #File = `file_id`
  )
## Info: applying file filter, keeping 67 of 71 files
## Info: applying file filter, keeping 67 of 67 files
## Info: aggregating vendor data table without units from 67 data file(s), including file info 'c(file_datetime, Analysis, `Identifier 1`, `Identifier 2`, `GC Method`, 
##     `Seed Oxidation`)'
data_table

Map peaks

### CHANGE MAPPING FILE NAME ###
metadata_samples <- read_excel(file.path("metadata", "20180310_GB_metadata.xlsx"), sheet = "files")

metadata_peak_maps <- read_excel(file.path("metadata","20180310_GB_metadata.xlsx"), sheet = "maps")

metadata_samples
data_table_with_peaks1 <- data_table %>%
    iso_add_metadata(metadata_samples, match_by = c(file_id)) %>%
    iso_map_peaks(metadata_peak_maps) 
## Info: metadata added to 1827 data rows, 0 left without metadata:
## - 64 metadata entries were mapped to 1827 data entires via column 'file_id'
## Info: 1422 peaks in 64 files were successfully assigned, 405 could not be assigned, and 506 were missing
# missing and unidentified peaks
data_table_with_peaks1 %>% filter(!is_identified | is_missing) %>% filter(process == "yes")
# confirmed peaks
data_table_with_peaks1 <- data_table_with_peaks1 %>%
    filter(is_identified, !is_missing, !is_ambiguous)

# missing and unidentified peaks
data_table_with_peaks <- data_table_with_peaks1 %>% filter(process == "yes")

data_table_with_peaks

Check stability of reference peaks

The second one is the one defined to be 0 permil (so will always be), the rest is relative to that peak.

p <- data_table_with_peaks %>% 
  filter(is_ref_peak == "yes") %>% 
  ggplot() + 
  aes(Nr., d13C, color = file_id) +
  geom_line() + 
  theme_bw()
ggplotly(p)

Add standard values

standards <- read_excel(file.path("metadata", "gc_irms_indiana_A6.xlsx"))
data_w_stds <- data_table_with_peaks %>% 
  filter(type == "standard", is_ref_peak == "no") %>% 
  left_join(standards, by = "compound") %>% 
  mutate(is_std = !is.na(true.d13C) | !is.na(true.d2H))

Process data

Focus on the analytes and calculate a few summary parameters we want to use later.

data_w_analyte_peaks <- 
  data_table_with_peaks %>% 
  # this is important so that the reference peaks are not caught up in the next set of calculations
  filter(is_ref_peak == "no") %>% 
  # for each analysis calculate averages across analysis
  group_by(file_id) %>% 
  mutate(
    ampl_sample_mean.mV = mean(`Ampl44`), ampl_sample_sd.mV = sd(`Ampl44`),
    area_sample_mean.Vs = mean(`Intensity 44`), area_sample_sd.Vs = sd(`Intensity 44`)
  )

Standards

standards <- read_excel(file.path("metadata", "gc_irms_indiana_A6.xlsx"))
kable(standards)
compound true.d2H true.d13C
C16 -9.1 -26.15
C17 -117.8 -31.88
C18 -52.0 -32.70
C19 -56.3 -31.99
C20 -89.7 -33.97
C21 -177.8 -28.83
C22 -81.3 -33.77
C23 -67.2 -33.37
C24 -29.7 -32.13
C25 -263.0 -28.46
C26 -45.9 -32.94
C27 -172.8 -30.49
C28 -36.8 -33.20
C29 -177.8 -29.10
C30 -213.6 -29.84 
data_w_stds1 <-
  data_w_analyte_peaks %>% 
  iso_add_standards(standards) 
## Info: added 21 standard entries to 564 out of 642 rows
data_w_stds <- data_w_stds1 %>% filter(type == "standard")

Visualize standards

v <- data_w_stds %>% 
  ggplot() +
  aes(x = true.d13C, y = d13C, color = file_id) + 
  geom_smooth(method = "lm", se = FALSE, alpha = 0.5) +
  geom_point() +
  theme_bw() +
  theme(legend.position = "none") 
ggplotly(v)
## We recommend that you use the dev version of ggplot2 with `ggplotly()`
## Install it with: `devtools::install_github('hadley/ggplot2')`

Calibration

data_w_calibs <- data_w_stds %>% 
  # prepare for calibration by defining the grouping column(s) and setting default parameters
  iso_prepare_for_calibration(group_by = c(file_id)) %>% 
  iso_set_default_process_parameters(delta_residual = resid_d13C) %>% 
  # run calibration
  iso_calibrate_delta(model = lm(`d13C` ~ true.d13C)) %>% 
  # pull out some columns we want generally available
  iso_unnest_calib_data(select = c(starts_with("Id"), ampl_sample_mean.mV)) 
## Info: preparing data for calibration by grouping based on 'file_id' and nesting the grouped datasets into 'calibration_data'
## Info: calculating delta calibration fits based on 1 model ('lm(d13C ~ true.d13C)') for 29 data group(s) in 'calibration_data' with filter 'is_standard'; storing residuals in 'resid_d13C.'
data_w_calibs %>% 
  iso_unnest_delta_calib_summary(keep_other_list_data = FALSE) %>% 
  kable(digits =3)
file_id Identifier 1 Identifier 2 ampl_sample_mean.mV r.squared adj.r.squared sigma statistic p.value df logLik AIC BIC deviance df.residual
590_A6_1_0uL.dxf A6 1.0uL 969.138 0.984 0.983 0.298 793.072 0 2 -2.060 10.120 12.244 1.156 13
591_A6_2_0uL.dxf A6 2.0uL 1955.459 0.996 0.995 0.154 3068.492 0 2 7.851 -9.702 -7.577 0.308 13
592_A6_1_0uL.dxf A6 1.0uL 985.106 0.987 0.986 0.264 995.667 0 2 -0.239 6.478 8.602 0.907 13
598_A6_0_5uL.dxf A6 0.5uL 394.483 0.916 0.908 0.634 109.424 0 2 -10.465 26.929 28.384 4.019 10
599_A6_1_0uL.dxf A6 1.0uL 1099.011 0.980 0.978 0.329 622.024 0 2 -3.519 13.038 15.162 1.404 13
600_A6_2_0uL.dxf A6 2.0uL 2268.953 0.993 0.992 0.198 1788.691 0 2 4.087 -2.175 -0.051 0.509 13
606_A6_0_5uL.dxf A6 0.5uL 396.842 0.921 0.913 0.637 116.625 0 2 -10.524 27.048 28.502 4.059 10
607_A6_1_5uL.dxf A6 1.5uL 1729.278 0.990 0.989 0.234 1265.227 0 2 1.602 2.796 4.920 0.709 13
608_A6_3_0uL.dxf A6 3.0uL 3132.502 0.999 0.998 0.087 9262.522 0 2 16.369 -26.738 -24.614 0.099 13
615_A6_0_5uL.dxf A6 0.5uL 455.763 0.949 0.945 0.511 243.367 0 2 -10.134 26.268 28.392 3.392 13
616_A6_1_5uL.dxf A6 1.5uL 693.561 0.959 0.954 0.404 209.650 0 2 -4.536 15.073 16.266 1.469 9
617_A6_3_0uL.dxf A6 3.0uL 1026.787 0.993 0.992 0.203 1429.118 0 2 3.210 -0.420 1.035 0.411 10
620_A6_0_5uL.dxf A6 0.5uL 136.079 0.824 0.802 1.299 37.387 0 2 -15.687 37.374 38.282 13.493 8
621_A6_1_5uL.dxf A6 1.5uL 755.950 0.952 0.947 0.507 216.766 0 2 -8.541 23.082 24.777 2.832 11
622_A6_3_0uL.dxf A6 3.0uL 1466.878 0.997 0.996 0.143 3151.419 0 2 7.952 -9.903 -8.208 0.224 11
629_A6_0_5uL.dxf A6 0.5uL 178.648 0.884 0.873 0.825 76.289 0 2 -13.626 33.252 34.706 6.807 10
630_A6_1_5uL.dxf A6 1.5uL 899.245 0.966 0.964 0.420 371.088 0 2 -7.191 20.382 22.506 2.291 13
631_A6_3_0uL.dxf A6 3.0uL 1540.907 0.993 0.992 0.198 1558.489 0 2 3.698 -1.396 0.299 0.431 11
640_A6_0_5uL.dxf A6 0.5uL 187.116 0.737 0.711 1.313 28.017 0 2 -19.200 44.401 45.856 17.238 10
641_A6_1_5uL.dxf A6 1.5uL 1179.021 0.978 0.977 0.337 586.638 0 2 -3.892 13.785 15.909 1.476 13
642_A6_3_0uL.dxf A6 3.0uL 1946.462 0.999 0.999 0.089 9589.068 0 2 16.054 -26.108 -23.984 0.103 13
644_A6_0_5uL.dxf A6 0.5uL 202.125 0.851 0.835 1.017 51.596 0 2 -14.688 35.377 36.570 9.306 9
645_A6_1_0uL.dxf A6 1.0uL 1048.358 0.977 0.975 0.349 471.183 0 2 -3.687 13.374 15.069 1.342 11
646_A6_2_0uL.dxf A6 2.0uL 2055.839 0.998 0.998 0.097 7928.598 0 2 14.776 -23.551 -21.427 0.122 13
647_A6_3_0uL.dxf A6 3.0uL 1691.292 0.998 0.998 0.097 7705.154 0 2 14.801 -23.603 -21.479 0.122 13
648_A6_0_5uL.dxf A6 0.5uL 156.361 0.942 0.935 0.503 145.239 0 2 -6.946 19.892 21.086 2.277 9
649_A6_1_0uL.dxf A6 1.0uL 1027.354 0.975 0.973 0.365 502.470 0 2 -5.099 16.197 18.322 1.733 13
650_A6_2_0uL.dxf A6 2.0uL 1776.514 0.998 0.998 0.105 6190.362 0 2 12.726 -19.453 -17.536 0.133 12
651_A6_3_0uL.dxf A6 3.0uL 1663.263 0.999 0.999 0.086 10332.483 0 2 16.583 -27.166 -25.042 0.096 13 
data_w_calibs %>% 
  iso_unnest_delta_calib_coefs(select = c(-statistic), keep_other_list_data = FALSE) %>% 
  # arrange by term and Analysis to get a quick idea of what the numbers across analyses
  arrange(term, file_id) %>% 
  kable(digits = 2)
file_id Identifier 1 Identifier 2 ampl_sample_mean.mV term estimate std.error p.value signif
590_A6_1_0uL.dxf A6 1.0uL 969.14 (Intercept) 9.11 1.07 0.00 ***
591_A6_2_0uL.dxf A6 2.0uL 1955.46 (Intercept) 9.79 0.56 0.00 ***
592_A6_1_0uL.dxf A6 1.0uL 985.11 (Intercept) 8.98 0.95 0.00 ***
598_A6_0_5uL.dxf A6 0.5uL 394.48 (Intercept) 11.59 3.20 0.00 **
599_A6_1_0uL.dxf A6 1.0uL 1099.01 (Intercept) 8.34 1.18 0.00 ***
600_A6_2_0uL.dxf A6 2.0uL 2268.95 (Intercept) 9.10 0.71 0.00 ***
606_A6_0_5uL.dxf A6 0.5uL 396.84 (Intercept) 13.03 3.21 0.00 **
607_A6_1_5uL.dxf A6 1.5uL 1729.28 (Intercept) 8.90 0.84 0.00 ***
608_A6_3_0uL.dxf A6 3.0uL 3132.50 (Intercept) 9.32 0.31 0.00 ***
615_A6_0_5uL.dxf A6 0.5uL 455.76 (Intercept) 7.26 1.84 0.00 **
616_A6_1_5uL.dxf A6 1.5uL 693.56 (Intercept) 11.18 2.25 0.00 ***
617_A6_3_0uL.dxf A6 3.0uL 1026.79 (Intercept) 10.32 0.83 0.00 ***
620_A6_0_5uL.dxf A6 0.5uL 136.08 (Intercept) 21.79 7.25 0.02 *
621_A6_1_5uL.dxf A6 1.5uL 755.95 (Intercept) 8.11 1.98 0.00 **
622_A6_3_0uL.dxf A6 3.0uL 1466.88 (Intercept) 10.60 0.56 0.00 ***
629_A6_0_5uL.dxf A6 0.5uL 178.65 (Intercept) 14.47 4.16 0.01 **
630_A6_1_5uL.dxf A6 1.5uL 899.25 (Intercept) 7.85 1.51 0.00 ***
631_A6_3_0uL.dxf A6 3.0uL 1540.91 (Intercept) 9.80 0.77 0.00 ***
640_A6_0_5uL.dxf A6 0.5uL 187.12 (Intercept) 12.59 6.62 0.09 .
641_A6_1_5uL.dxf A6 1.5uL 1179.02 (Intercept) 8.32 1.21 0.00 ***
642_A6_3_0uL.dxf A6 3.0uL 1946.46 (Intercept) 10.58 0.32 0.00 ***
644_A6_0_5uL.dxf A6 0.5uL 202.12 (Intercept) 18.79 5.67 0.01 **
645_A6_1_0uL.dxf A6 1.0uL 1048.36 (Intercept) 8.71 1.36 0.00 ***
646_A6_2_0uL.dxf A6 2.0uL 2055.84 (Intercept) 10.15 0.35 0.00 ***
647_A6_3_0uL.dxf A6 3.0uL 1691.29 (Intercept) 9.71 0.35 0.00 ***
648_A6_0_5uL.dxf A6 0.5uL 156.36 (Intercept) 12.04 2.80 0.00 **
649_A6_1_0uL.dxf A6 1.0uL 1027.35 (Intercept) 8.40 1.32 0.00 ***
650_A6_2_0uL.dxf A6 2.0uL 1776.51 (Intercept) 10.00 0.39 0.00 ***
651_A6_3_0uL.dxf A6 3.0uL 1663.26 (Intercept) 10.72 0.31 0.00 ***
590_A6_1_0uL.dxf A6 1.0uL 969.14 true.d13C 0.97 0.03 0.00 ***
591_A6_2_0uL.dxf A6 2.0uL 1955.46 true.d13C 0.98 0.02 0.00 ***
592_A6_1_0uL.dxf A6 1.0uL 985.11 true.d13C 0.96 0.03 0.00 ***
598_A6_0_5uL.dxf A6 0.5uL 394.48 true.d13C 1.07 0.10 0.00 ***
599_A6_1_0uL.dxf A6 1.0uL 1099.01 true.d13C 0.94 0.04 0.00 ***
600_A6_2_0uL.dxf A6 2.0uL 2268.95 true.d13C 0.96 0.02 0.00 ***
606_A6_0_5uL.dxf A6 0.5uL 396.84 true.d13C 1.11 0.10 0.00 ***
607_A6_1_5uL.dxf A6 1.5uL 1729.28 true.d13C 0.96 0.03 0.00 ***
608_A6_3_0uL.dxf A6 3.0uL 3132.50 true.d13C 0.97 0.01 0.00 ***
615_A6_0_5uL.dxf A6 0.5uL 455.76 true.d13C 0.92 0.06 0.00 ***
616_A6_1_5uL.dxf A6 1.5uL 693.56 true.d13C 1.03 0.07 0.00 ***
617_A6_3_0uL.dxf A6 3.0uL 1026.79 true.d13C 1.00 0.03 0.00 ***
620_A6_0_5uL.dxf A6 0.5uL 136.08 true.d13C 1.41 0.23 0.00 ***
621_A6_1_5uL.dxf A6 1.5uL 755.95 true.d13C 0.94 0.06 0.00 ***
622_A6_3_0uL.dxf A6 3.0uL 1466.88 true.d13C 1.01 0.02 0.00 ***
629_A6_0_5uL.dxf A6 0.5uL 178.65 true.d13C 1.16 0.13 0.00 ***
630_A6_1_5uL.dxf A6 1.5uL 899.25 true.d13C 0.93 0.05 0.00 ***
631_A6_3_0uL.dxf A6 3.0uL 1540.91 true.d13C 0.98 0.02 0.00 ***
640_A6_0_5uL.dxf A6 0.5uL 187.12 true.d13C 1.12 0.21 0.00 ***
641_A6_1_5uL.dxf A6 1.5uL 1179.02 true.d13C 0.94 0.04 0.00 ***
642_A6_3_0uL.dxf A6 3.0uL 1946.46 true.d13C 1.00 0.01 0.00 ***
644_A6_0_5uL.dxf A6 0.5uL 202.12 true.d13C 1.29 0.18 0.00 ***
645_A6_1_0uL.dxf A6 1.0uL 1048.36 true.d13C 0.95 0.04 0.00 ***
646_A6_2_0uL.dxf A6 2.0uL 2055.84 true.d13C 0.99 0.01 0.00 ***
647_A6_3_0uL.dxf A6 3.0uL 1691.29 true.d13C 0.98 0.01 0.00 ***
648_A6_0_5uL.dxf A6 0.5uL 156.36 true.d13C 1.07 0.09 0.00 ***
649_A6_1_0uL.dxf A6 1.0uL 1027.35 true.d13C 0.94 0.04 0.00 ***
650_A6_2_0uL.dxf A6 2.0uL 1776.51 true.d13C 0.98 0.01 0.00 ***
651_A6_3_0uL.dxf A6 3.0uL 1663.26 true.d13C 1.01 0.01 0.00 ***

Parameters

data_params <- data_w_calibs %>% 
  # pull out remaining columns we want available (some might already be pulled out but that's okay)
  #Note: Garrett changed "Preparation" to "Seed Oxidation" - we may want a different shape variable
  iso_unnest_calib_data(select = c(file_datetime, ampl_sample_mean.mV, `Seed Oxidation`, is_standard)) 
  
# visualize the delta calibration fits
data_params %>% 
  #NEED TO FILTER OUT SAMPLES iso_filter_files(`Identifier 1` == "A5") 
  iso_visualize_delta_calib_fits(x = file_id, color = `Identifier 2`, size = ampl_sample_mean.mV,
                                 include_from_summary = c(adj.r.squared, deviance)) + labs(title = "parameters vs. analysis")

data_params %>% 
  iso_visualize_delta_calib_fits(x = ampl_sample_mean.mV, color = `Identifier 2`, size = ampl_sample_mean.mV,
                                 include_from_summary = c(adj.r.squared, deviance)) + labs(title = "parameters vs. amplitude") 

data_params %>% 
  iso_visualize_delta_calib_fits(x = file_datetime, color = `Identifier 2`, size = ampl_sample_mean.mV,
                                 include_from_summary = c(adj.r.squared, deviance)) + labs(title = "parameters vs. time")

# turn the last plot into an interactive one
ggplotly(ggplot2::last_plot() + theme(legend.position = "none"))
## We recommend that you use the dev version of ggplot2 with `ggplotly()`
## Install it with: `devtools::install_github('hadley/ggplot2')`

Compounds

data_w_calibs %>% 
  # pull out relevant data
  iso_unnest_calib_data(select = everything()) %>% 
  filter(is_standard) %>% 
  # calculate deviation from means
  group_by(file_id) %>% 
  mutate(
    `Var: residual d13C [permil]` = resid_d13C,
    `Var: area diff from mean [%]` = (`Intensity 44`/mean(`Intensity 44`) - 1) * 100,
    `Var: amplitude diff from mean [%]` = (`Ampl44`/mean(`Ampl44`) - 1) * 100
  ) %>%
  # visualize
  iso_visualize_data(x = compound, y = starts_with("Var"), group = file_id, color = `Identifier 2`)

ggplotly(ggplot2::last_plot())
## We recommend that you use the dev version of ggplot2 with `ggplotly()`
## Install it with: `devtools::install_github('hadley/ggplot2')`

Apply correction to Adducts